function Gamma_ = MPC_2Huggett(C_, S_, T, N, incomeStruct, Xstruct)
%{
%%-------------------------------------------------------------------------
Solves the Feynman-Kac MPC formula for the 2-state Huggett economy
%%-------------------------------------------------------------------------
%}    
    
%%-------------------------------------------------------------------------
%Unpack Structs / Setup
    
    lambda_L = incomeStruct.lambda_L;
    lambda_H = incomeStruct.lambda_H;
    
    I = Xstruct.I;
    xjump = Xstruct.xjump;
    
    tjump = T / N;
    Gamma_ = zeros(I, 2, N+1); 
    
%%-------------------------------------------------------------------------


%%-------------------------------------------------------------------------
%Solve PDE for cumulative consumption 

    %Initialize
    Aswitch = [-speye(I)*lambda_L, speye(I)*lambda_L; speye(I)*lambda_H, -speye(I)*lambda_H];
    Ib = (S_ < 0);
    If = (S_ > 0);
    I0 = (S_ == 0);
    
    X = -Ib.*S_ / xjump;
    Y = -If.*S_ / xjump + Ib.*S_ / xjump;
    Z = If.*S_ / xjump;

    A1=spdiags(Y(:,1),0,I,I)+spdiags(X(2:I,1),-1,I,I)+spdiags([0;Z(1:I-1,1)],1,I,I);
    A2=spdiags(Y(:,2),0,I,I)+spdiags(X(2:I,2),-1,I,I)+spdiags([0;Z(1:I-1,2)],1,I,I);
    A = [A1,sparse(I,I);sparse(I,I),A2] + Aswitch;

    B = (1/tjump)*speye(2*I) - A;

    C_stacked = [C_(:,1); C_(:,2)];
    
    %Calculate cumulative consumption
    for n = N+1:-1:2
        Gamma_stacked = [Gamma_(:,1, n); Gamma_(:, 2, n)];

        b = C_stacked + Gamma_stacked/tjump;
        Gamma_stacked = B\b;

        Gamma_(:, 1, n-1) = Gamma_stacked(1:I);
        Gamma_(:, 2, n-1) = Gamma_stacked(I+1:2*I);       
    end    
%%-------------------------------------------------------------------------
    

end